from pathlib import Path
import subprocess
from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u
import os
from astropy.table import Table, vstack, hstack
import numpy as np
from astropy.io import fits
from imp import reload
import matplotlib.pyplot as plt
import os
from astropy.coordinates import SkyCoord
import astropy.units as u
import pandas as pd
from astropy.wcs import WCS
from astropy.stats import sigma_clipped_stats
from astropy.stats import sigma_clip
from astropy.time import Time
import seaborn as sns
from datetime import datetime, timedelta
import matplotlib.colors as colors
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.ticker as ticker
import sep
from loguru import logger
import warnings

from tools.util_statck import run_command
warnings.filterwarnings('ignore')
def read_image_list(file_path):
    """
    读取指定文件并遍历每一行数据。

    参数:
    file_path (str): 包含图像路径的文件路径。

    返回:
    list: 有效的图像路径列表。
    """
    image_paths = []
    if os.path.exists(file_path)==False:
        logger.info(f"file_path:{file_path} not exists")
        return []
    try:
        with open(file_path, 'r') as file:
            for line in file:
                # 去除行末的换行符和空格
                image_path = line.strip()
                if image_path:  # 确保不是空行
                    image_paths.append(image_path)  # 添加到列表中
                    # 在这里可以对每个 image_path 进行处理
                    logger.info(image_path)  # 示例：打印每个图像路径
    except FileNotFoundError:
        logger.info(f"文件未找到: {file_path}")
    except Exception as e:
        logger.info(f"发生错误: {e}")

    return image_paths

def save_region(target_ra, target_dec, name='', filename="region.reg"):
    """
    保存目标 RA 和 Dec 为区域文件格式.

    参数:
    target_ra (float): 目标赤经
    target_dec (float): 目标赤纬
    name (str): 目标名称
    filename (str): 保存的文件名，默认为 "region.reg"
    """
    # 创建区域字符串
    region_string = f"fk5; circle({target_ra}, {target_dec}, 10\") # {name}\n"

    # 将区域字符串写入文件
    with open(filename, "w") as file:
        file.write(region_string)

    logger.info(f"区域已保存到 {filename}")

    
def TargetXCatalogue(ra, dec, catalog, radius=1.2, RaKeyword='ALPHA_J2000', DecKeyword='DELTA_J2000'):
    RaDecTarget = SkyCoord(ra=ra*u.degree, dec=dec*u.degree)
    RaList = catalog[RaKeyword]
    DecList = catalog[DecKeyword]
    
    try:
        RaDecList = SkyCoord(ra=RaList, dec=DecList)
    except Exception as e:
        RaDecList = SkyCoord(ra=RaList * u.deg, dec=DecList * u.deg)
        
    idx, d2d, d3d = RaDecTarget.match_to_catalog_3d(RaDecList, nthneighbor=1)
    
    MaxSep = radius * u.arcsec
    sep_constraint = d2d < MaxSep
    
    if sep_constraint[0]:  # 如果有匹配
        return catalog[idx]  # 返回匹配的星表
    else:
        return None  # 返回空
def convertHMSDMS(ra,dec):
    coor=SkyCoord(ra*u.deg,dec*u.deg)
    rahms='%02d:%02d:%05.2f'%(coor.ra.hms.h,coor.ra.hms.m,coor.ra.hms.s)
    DecDMSD = coor.dec.dms.d
    if DecDMSD == 0:
        if '-' in str(DecDMSD):
            decdms='-00:%02d:%02d'%(abs(coor.dec.dms.m),abs(coor.dec.dms.s))
        else:
            decdms='%02d:%02d:%02d'%(coor.dec.dms.d,abs(coor.dec.dms.m),abs(coor.dec.dms.s))
    elif -10 < DecDMSD < 0:
        decdms='%03d:%02d:%02d'%(coor.dec.dms.d,abs(coor.dec.dms.m),abs(coor.dec.dms.s))
    else:
        decdms='%02d:%02d:%02d'%(coor.dec.dms.d,abs(coor.dec.dms.m),abs(coor.dec.dms.s))
    #logger.info(rahms,decdms)
    return rahms,decdms

def replace_suffix_pathlib(cata_path, mag_type):
    p = Path(cata_path)
    new_name = f"{p.stem}_{mag_type}_SourceDetection.reg"
    return p.with_name(new_name)

def CatalogueRegion(CataPath, MagType='AUTO'):
    """
    Generate a DS9 region file from a FITS catalog.

    Parameters
    ----------
    CataPath : str
        Path to the input catalog (FITS format).
    MagType : str, optional
        Used to name the output region file. Default is 'AUTO'.
    """

    # Define output region file name
    RagionPath = replace_suffix_pathlib(CataPath,MagType) #CataPath.replace('.cat', f'_{MagType}_SourceDetection.reg')
    logger.info(f"Writing region file: {RagionPath}")

    # Open the FITS file using a context manager
    with fits.open(CataPath) as hdul:
        # Attempt to read from extension [2]. If it fails, fallback to [1].
        try:
            cata_data = Table(hdul[2].data)
        except Exception:
            cata_data = Table(hdul[1].data)

    # Grab RA/Dec columns
    RaList = cata_data['ALPHA_J2000']
    DecList = cata_data['DELTA_J2000']
    snrList = cata_data['SNR_WIN']
    mag_AUTOList = cata_data['MAG_AUTO']

    # Write the region file only once
    with open(RagionPath, 'w') as reglist:
        reglist.write('# Region file for DS9\n')
        reglist.write('global color=red font="helvetica 10 normal" '
                      'select=1 edit=1 move=1 delete=1 include=1 '
                      'fixed=0 source=1\n')
        reglist.write('fk5\n')
        for ra_val, dec_val,snr_val,mag_auto in zip(RaList, DecList,snrList,mag_AUTOList):
            region_line = f'circle({ra_val},{dec_val},1")  '
            region_line += '# color=red font="helvetica 12 normal roman"'
            region_line += f"text={{SNR={snr_val:.2f} MAG_AUTO={mag_auto:.2f}}}"
            reglist.write(region_line + "\n")


def remove_flxscale_from_fits(fits_file):
    """
    从 FITS 文件中删除 FLXSCALE 关键字并覆盖原文件.

    参数:
    fits_file (str): 输入的 FITS 文件路径
    """
    logger.info(f"从 FITS:{fits_file} 文件中删除 FLXSCALE 关键字并覆盖原文件.")
    # 打开 FITS 文件
    with fits.open(fits_file, mode='update') as hdul:
        # 遍历每个头部
        for hdu in hdul:
            # 检查并删除 FLXSCALE 关键字
            if 'FLXSCALE' in hdu.header:
                logger.info(f"Removing FLXSCALE from header of HDU {hdu.name}")
                del hdu.header['FLXSCALE']
        
        # 保存修改后的文件，覆盖原文件
        hdul.flush()  # 确保更改被写入
        logger.info(f"Modified FITS file saved as {fits_file}")
        
def imageCenter(wcsImg):
    """
    find the celestical center of a given image wcs
    """
    ximg, yimg = wcsImg.pixel_shape
    ximgCen, yimgCen = 0.5*(ximg+1), 0.5*(yimg+1)
    raCen, decCen    = wcsImg.all_pix2world(ximgCen, yimgCen,1)
    raCenS, decCenS  = convertHMSDMS(raCen, decCen)
    return raCenS, decCenS



def WcsAlignFunc(imgList, WcsAlignedImgPathFolder, WcsSuffix="_WcsAligned.fits", SwarpConfig='./config/default_resampling.swarp', COMBINE="N", RESAMPLE="Y"):
    
    os.makedirs(WcsAlignedImgPathFolder, exist_ok=True)
    # 定义一个函数 WcsAlignFunc，接受多个参数，包括图像列表、输出文件夹、文件后缀、Swarp 配置文件路径、合并选项和重采样选项。

    SwarpWcsAlign1 = "swarp %s -c %s -CENTER %s,%s -PIXEL_SCALE %f -IMAGE_SIZE %d,%d -SUBTRACT_BACK N -VERBOSE_TYPE FULL "
    # 定义 Swarp 命令的第一部分，包含输入图像、配置文件、中心坐标、像素比例、图像大小等参数。

    SwarpWcsAlign2 = "-IMAGEOUT_NAME %s -WEIGHTOUT_NAME %s -RESAMPLE_DIR %s -COMBINE %s -RESAMPLE %s -FSCALE_KEYWORD FLXSCALE_TEMP"
    # 定义 Swarp 命令的第二部分，包含输出图像名称、权重图像名称、重采样目录、合并和重采样选项及临时缩放关键字。

    SwarpWcsAlign = SwarpWcsAlign1 + SwarpWcsAlign2
    # 将两部分命令合并为一个完整的 Swarp 命令字符串。

    SwarpConfig = SwarpConfig
    # 将传入的 Swarp 配置文件路径赋值给 SwarpConfig 变量（此行实际上没有改变内容）。

    imgoutList = []
    # 初始化一个空列表，用于存储输出图像的路径。

    pixScale, imgSize, cenImg = [], [], []
    # 初始化三个空列表，分别用于存储像素比例、图像大小和图像中心坐标。

    for i in range(2):
        # 循环处理前两个图像（假设 imgList 至少有两个图像）。

        iimg = imgList[i]
        # 获取当前图像的文件名。

        iimgout = iimg.split("/")[-1][:-5] + WcsSuffix
        # 从当前图像的文件名中提取文件名（去掉路径和文件扩展名），并添加后缀。

        iimgout = os.path.join(WcsAlignedImgPathFolder, iimgout)
        # 将输出图像的文件名与输出文件夹路径结合，得到完整的输出路径。

        imgoutList += [iimgout]
        # 将输出图像的路径添加到 imgoutList 列表中。

        iwcs = WCS(iimg)
        # 使用 WCS 库加载当前图像的 WCS 信息。

        header = fits.getheader(iimg)
        # 获取当前图像的 FITS 头信息。

        band = header['FILTER']
        # 从头信息中提取滤光片信息（band）。

        ipixScale = 3600 * np.mean(np.sqrt(np.sum(iwcs.pixel_scale_matrix**2, axis=1)))
        # 计算图像的像素比例，单位为弧秒。

        if band == 'u' or band == 'v' or band == 'g' or band == 'r':
            # 如果滤光片是 u、v、g 或 r。

            ipixScale = 0.429
            # 设置像素比例为 0.429。

            pixScale += [float("%.3f" % ipixScale)]
            # 将像素比例添加到 pixScale 列表中，并格式化为三位小数。

        elif band == 'i' or band == 'z':
            # 如果滤光片是 i 或 z。

            ipixScale = 0.286
            # 设置像素比例为 0.286。

            pixScale += [float("%.3f" % ipixScale)]
            # 将像素比例添加到 pixScale 列表中，并格式化为三位小数。

        # logger.info('Pixel scale is %s' % (pixScale))
        # （注释掉的日志记录）记录当前像素比例的信息。

        imgSize += [iwcs.pixel_shape]
        # 将当前图像的像素形状（大小）添加到 imgSize 列表中。

        cenImg += [imageCenter(iwcs)]
        # 计算当前图像的中心坐标，并添加到 cenImg 列表中。

    for i in range(0, 2):
        # 循环处理前两个图像。

        iimg = imgList[i]
        # 获取当前图像的文件名。

        resam_dir, _ = os.path.split(iimgout)
        # 从输出图像路径中提取重采样目录。

        iimgout = imgoutList[i]
        # 获取当前图像的输出路径。

        wghTMP = iimgout[:-5] + "_wgh.fits"
        # 生成权重图像的临时文件名。

        comdSwarpIn = SwarpWcsAlign % (iimg, SwarpConfig, cenImg[1][0], cenImg[1][1], pixScale[i],
                                        imgSize[1][0], imgSize[1][1], iimgout, wghTMP, resam_dir, COMBINE, RESAMPLE)
        # 格式化完整的 Swarp 命令，插入各个参数。

        logger.info(comdSwarpIn)
        # 记录生成的 Swarp 命令。

        os.system(comdSwarpIn)
        # 使用系统命令执行 Swarp 命令。

    return imgoutList
    # 返回输出图像的路径列表。



def psfex_model(catalog_name,psfex_config,snr_low=10, kernel_size_factor=7):
    """
    创建PSFEx模型，并返回模型文件路径。
    """
    psfname = Path(catalog_name).with_suffix('').name  # 从星表文件名派生PSF文件名
    psf_dir = Path(catalog_name).parent / psfname  # PSF模型存储目录

    # 如果PSF模型目录不存在，则创建它
    if not psf_dir.exists():
        psf_dir.mkdir(exist_ok=True)  # 
        #创建PSF模型目录
        
    cata_data = Table.read(catalog_name, hdu=2)  # 读取星表数据
    # 选择信噪比大于snr_low且标志位为0的对象
    original_catalog_select = cata_data[(cata_data['SNR_WIN'] > snr_low) & (cata_data['FLAGS'] == 0)]
    # 计算选定对象的FWHM的中位数
    median_fwhm = sigma_clipped_stats(original_catalog_select['FWHM_IMAGE'], sigma=3)[1]
    kernel_size = int(median_fwhm) * kernel_size_factor  # 计算内核大小

    checkimage_dir = psf_dir

    
    # 构建PSFEx命令
    psfex_command = (
        f'psfex {catalog_name} '
        f'-c {psfex_config} '
        f'-PSF_SIZE {kernel_size},{kernel_size} '
        f'-PSF_DIR {psf_dir} '
        f'-OUTCAT_NAME {psf_dir}/psfex_out.ldac '
        f'-XML_NAME {psf_dir}/psfex.xml '
        f"-CHECKPLOT_TYPE FWHM,ELLIPTICITY,COUNTS,CHI2,RESIDUALS "
        f"-CHECKPLOT_NAME {checkimage_dir}/fwhm,{checkimage_dir}/ellipticity,{checkimage_dir}/counts,{checkimage_dir}/chi2,{checkimage_dir}/resi "
        f"-CHECKIMAGE_TYPE CHI,PROTOTYPES,SAMPLES,RESIDUALS,SNAPSHOTS "
        f"-CHECKIMAGE_NAME {','.join([os.path.join(checkimage_dir, name) for name in ['chi.fits', 'proto.fits', 'samp.fits', 'resi.fits', 'snap.fits']])} "
    )
    logger.info(psfex_command)
    run_command(psfex_command)  # 执行PSFEx命令
    psf_model_path=psf_dir / f"{psfname}.psf"
    if os.path.exists(psf_model_path)==False:
        raise
    return psf_model_path  # 返回PSF模型文件路径
    
    
def SexPhotoTwoImg1m6(SciimgPath, SourceDetectionPath, Psf_path, SexConfig=None, SexParam=None, CataSuffix=None, PsfPhot=False,DETECT_MINAREA=5,DETECT_THRESH=1.5,ANALYSIS_THRESH=1.5):
    """
    使用 SExtractor 对两个图像进行光度测量。

    参数:
    SciimgPath (str): 科学图像的路径
    SourceDetectionPath (str): 源检测图像的路径
    PsfSuffix (str): PSF 文件的后缀
    SexConfig (str, optional): SExtractor 配置文件路径。默认为 None。
    SexParam (str, optional): SExtractor 参数文件路径。默认为 None。
    CataSuffix (str, optional): 输出目录文件的后缀。默认为 None。
    PsfPhot (bool, optional): 是否进行 PSF 光度测量。默认为 False。

    返回:
    str: 输出目录文件的路径
    """
    
    if PsfPhot == False:
        # 如果不进行 PSF 光度测量
        sciimgpath1 = SourceDetectionPath
        sciimgpath2 = SciimgPath
        # logger.info(sciimgpath1, sciimgpath2)
        psfpath = None  # PSF 文件路径设为 None
    elif PsfPhot == True:
        # 如果进行 PSF 光度测量
        sciimgpath1 = SourceDetectionPath
        sciimgpath2 = SciimgPath
        psfpath = Psf_path  # 生成 PSF 文件路径

    # 读取科学图像的头信息
    sciimghdu = fits.open(sciimgpath2)
    header = sciimghdu[0].header
    gain = header['GAIN']  # 获取增益值
    saturate = header['SATURATE']  # 获取饱和水平
    FilenameInfo = sciimgpath2.split('.')[0]  # 获取文件名信息
    CatalogueName = FilenameInfo + '%s.cat' % CataSuffix  # 生成目录文件名

    try:
        # 尝试从头信息中获取视宁度
        seeing = float(header['SEEING'])
    except:
        # 如果头信息中没有视宁度，使用默认值 1.2
        seeing = 1.2

    # 构建 SExtractor 命令
    sexComd1 = "sex %s %s -c %s -PARAMETERS_NAME %s -CATALOG_NAME %s" % (sciimgpath1, sciimgpath2, SexConfig, SexParam, CatalogueName)
    sexComd2 = f" -DETECT_MINAREA {DETECT_MINAREA} -DETECT_THRESH {DETECT_THRESH} -ANALYSIS_THRESH {ANALYSIS_THRESH}"
    sexComd3 = " -SATUR_LEVEL %s -SEEING_FWHM %s -GAIN %s -CHECKIMAGE_TYPE NONE" % (saturate, seeing, gain)
    sexComd4 = " -PSF_NAME %s" % psfpath

    if psfpath is None:
        # 如果没有 PSF 文件路径，构建不包含 PSF 的命令
        sexComd = sexComd1 + sexComd2 + sexComd3
        logger.info(sexComd)  # 记录命令
        os.system(sexComd)  # 执行命令
        logger.info('Photometry of %s is OK!' % sciimgpath2)  # 记录光度测量完成
        return CatalogueName  # 返回目录文件名
    elif psfpath is not None:
        # 如果有 PSF 文件路径，构建包含 PSF 的命令
        sexComdPsf = sexComd1 + sexComd2 + sexComd3 + sexComd4
        logger.info(sexComdPsf)
        os.system(sexComdPsf)  # 执行命令
        logger.info('Photometry of %s is OK!' % sciimgpath2)  # 记录光度测量完成
        return CatalogueName  # 返回目录文件名
    
import math

def convert_to_unit_time_magnitude(m, t, sigma_m=0.0, sigma_t=0.0):
    """
    将总星等转换为单位时间星等，并计算相应的误差。

    参数:
    m (float): 总星等
    t (float): 曝光时间（秒）
    sigma_m (float): 总星等的误差（默认为0.0）
    sigma_t (float): 曝光时间的误差（默认为0.0）

    返回:
    tuple: (单位时间星等, 单位时间星等的误差)
    """
    # 计算单位时间星等
    m_per_sec = m + 2.5 * math.log10(t)
    
    # 计算单位时间星等的误差
    if sigma_t == 0.0:
        sigma_m_per_sec = sigma_m
    else:
        term1 = (2.5 * sigma_t) / (t * math.log(10))
        sigma_m_per_sec = math.sqrt(sigma_m**2 + term1**2)
    
    return m_per_sec, sigma_m_per_sec

# if __name__ == "__main__":
#     # 示例使用
#     m = -6.57  # 总星等
#     t = 300.0  # 曝光时间（秒）
#     sigma_m = 0.01  # 总星等的误差
#     sigma_t = 0.1  # 曝光时间的误差

#     m_per_sec, sigma_m_per_sec = convert_to_unit_time_magnitude(m, t, sigma_m, sigma_t)

#     logger.info(f"单位时间星等: {m_per_sec:.4f}")
#     logger.info(f"单位时间星等的误差: {sigma_m_per_sec:.4f}")




def read_header_keys(fits_file, keys):
    """
    读取 FITS 文件头信息中的指定关键字。

    参数:
    fits_file (str): FITS 文件的路径
    keys (list of str): 要读取的关键字列表

    返回:
    dict: 包含关键字及其对应值的字典
    """
    header_values = {}
    
    # 读取 FITS 文件
    with fits.open(fits_file) as hdul:
        # 获取主头信息（通常是第一个 HDU）
        header = hdul[0].header
        
        # 遍历关键字列表，读取对应的值
        for key in keys:
            value = header.get(key, None)
            header_values[key] = value
    
    return header_values

def ensure_directories(*dirs):
    """
    确保所有指定的目录存在。
    """
    for directory in dirs:
        Path(directory).mkdir(parents=True, exist_ok=True)
        
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

def plot_calmag_with_error_bars(file_path,png_):
    """
    读取 CSV 文件并绘制 CALMAG 和 Sigma_MAG_PER_SEC 的图形。

    参数:
    file_path (str): CSV 文件的路径。
    """
    # 读取 CSV 文件
    data = pd.read_csv(file_path)

    # 提取时间和相关数据
    data['DATE'] = pd.to_datetime(data['DATE'])  # 将 DATE 列转换为 datetime 格式
    calmag = data['CALMAG']
    sigma_mag_per_sec = data['Sigma_MAG_PER_SEC']

    # 绘图
    plt.figure(figsize=(10, 6))
    plt.errorbar(data['DATE'], calmag, yerr=sigma_mag_per_sec, fmt='o', capsize=5, label='CALMAG with Error Bars', color='blue')
    
    # 添加文本标签
    for i in range(len(data)):
        plt.text(data['DATE'].iloc[i], calmag.iloc[i], f"{calmag.iloc[i]:.2f} ± {sigma_mag_per_sec.iloc[i]:.2f}", 
                 fontsize=12, ha='center', va='bottom')

    # 图形美化
    plt.title('CALMAG vs Time with Error Bars')
    plt.xlabel('Time')
    plt.ylabel('CALMAG')

    # 格式化 x 轴，保留时分秒
    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
    plt.gca().xaxis.set_major_locator(mdates.MinuteLocator(interval=1))  # 根据数据频率调整间隔
    plt.gcf().autofmt_xdate()  # 自动旋转日期标签

    
    plt.legend()
    plt.grid()
    plt.tight_layout()
    
    # 翻转 y 轴
    plt.gca().invert_yaxis()

    # 显示图形
    
    # plt.show()
    
    plt.savefig(png_,dpi=300)

# 使用示例
# plot_calmag_with_error_bars('your_file_path.csv')
import math
from typing import Tuple, Dict, Optional, Any
class MatrixPartitioner:
    """
    一个用于划分大矩阵为子矩阵并根据坐标查找所在子矩阵的工具类。
    按照列优先、从下到上的顺序编号子矩阵。
    """

    def __init__(self, row: int, col: int):
        """
        初始化 MatrixPartitioner 类。

        参数:
        - row (int): 大矩阵的行数。
        - col (int): 大矩阵的列数。
        """
        self.row = row
        self.col = col

        # 根据矩阵尺寸自动设置子矩阵划分方式
        if self.row ==  6160 and self.col ==6144:
            self.subrow = 2
            self.subcol = 2
        elif self.row == 9232 and self.col == 9216:
            self.subrow = 2
            self.subcol = 8
        else:
            raise ValueError(f"不支持的矩阵尺寸: {self.row}x{self.col}")

        self.rows_per_submatrix = row // self.subrow
        self.cols_per_submatrix = col // self.subcol

        self._validate_parameters()
        self.submatrices = self._generate_submatrix_offsets()

    def _validate_parameters(self):
        """
        验证输入参数的有效性。
        """
        if self.subrow <= 0 or self.subcol <= 0:
            raise ValueError("subrow 和 subcol 必须是正整数。")
        if self.row <= 0 or self.col <= 0:
            raise ValueError("row 和 col 必须是正整数。")

    def _generate_submatrix_offsets(self) -> Dict[int, Dict[str, Any]]:
        """
        生成子矩阵的偏移量和尺寸信息。
        按照列优先、从下到上的顺序编号子矩阵。

        返回:
        - dict: 一个字典，键为子矩阵的编号（从1开始），值为包含偏移量和尺寸的字典。
        """
        submatrices = {}
        index = 1
        for j in range(self.subcol):  # 水平方向（列）优先
            for i in range(self.subrow):  # 垂直方向（行）从下到上
                y_start = i * self.rows_per_submatrix
                x_start = j * self.cols_per_submatrix

                # 处理最后一行或列的子矩阵，确保包含所有剩余部分
                if i == self.subrow - 1:
                    y_end = self.row
                else:
                    y_end = y_start + self.rows_per_submatrix

                if j == self.subcol - 1:
                    x_end = self.col
                else:
                    x_end = x_start + self.cols_per_submatrix

                submatrices[index] = {
                    'x_offset': x_start,
                    'y_offset': y_start,
                    'width': x_end - x_start,
                    'height': y_end - y_start,
                    'x_end': x_end,
                    'y_end': y_end
                }
                index += 1
        return submatrices

    def get_submatrix_offsets(self) -> Dict[int, Dict[str, Any]]:
        """
        获取所有子矩阵的偏移量和尺寸信息。

        返回:
        - dict: 子矩阵偏移量和尺寸信息。
        """
        return self.submatrices

    def find_submatrix(self, x: int, y: int) -> int:
        """
        根据给定的 (x, y) 坐标，判断其位于哪个子矩阵中。

        参数:
        - x (int): 水平方向的坐标（列索引，从0开始）。
        - y (int): 垂直方向的坐标（行索引，从0开始）。

        返回:
        - int: 子矩阵的编号（从1开始）。
        """
        if not (0 <= x < self.col) or not (0 <= y < self.row):
            raise ValueError("坐标 (x, y) 超出矩阵范围。")

        submatrix_row = y // self.rows_per_submatrix
        submatrix_col = x // self.cols_per_submatrix

        # 处理边界情况
        if submatrix_row >= self.subrow:
            submatrix_row = self.subrow - 1
        if submatrix_col >= self.subcol:
            submatrix_col = self.subcol - 1

        # 按照列优先、从下到上的顺序编号
        submatrix_index = submatrix_col * self.subrow + submatrix_row + 1
        return int(submatrix_index)

    def find_submatrix_with_range(self, x: int, y: int) -> Dict[str, Any]:
        """
        根据给定的 (x, y) 坐标，判断其位于哪个子矩阵中，并返回子矩阵的编号及其范围。

        参数:
        - x (int): 水平方向的坐标（列索引，从0开始）。
        - y (int): 垂直方向的坐标（行索引，从0开始）。

        返回:
        - dict: 包含子矩阵编号和范围的信息。
        """
        try:
            submatrix_index = self.find_submatrix(x, y)

            submatrix_info = self.submatrices[submatrix_index]
            return {
                'submatrix_index': submatrix_index,
                'y_range': (submatrix_info['y_offset'], submatrix_info['y_end'] - 1),
                'x_range': (submatrix_info['x_offset'], submatrix_info['x_end'] - 1)
            }
        except:
            return None

    def extract_submatrix_data(self, matrix: Any, submatrix_index: int) -> Any:
        """
        根据子矩阵编号，从大矩阵中提取子矩阵的数据。

        参数:
        - matrix (Any): 大矩阵的数据，通常是 NumPy 数组或类似的数据结构。
        - submatrix_index (int): 子矩阵的编号。

        返回:
        - Any: 子矩阵的数据。
        """
        if submatrix_index not in self.submatrices:
            raise ValueError(f"子矩阵编号 {submatrix_index} 不存在。")

        submatrix_info = self.submatrices[submatrix_index]
        y_start = submatrix_info['y_offset']
        y_end = submatrix_info['y_end']
        x_start = submatrix_info['x_offset']
        x_end = submatrix_info['x_end']

        return matrix[y_start:y_end, x_start:x_end]

    def __str__(self):
        """
        返回子矩阵的详细信息。
        """
        info = f"MatrixPartitioner: {self.row}x{self.col} split into {self.subrow}x{self.subcol} submatrices.\n"
        for idx, sub in self.submatrices.items():
            info += f" Submatrix {idx}: x_offset={sub['x_offset']}, y_offset={sub['y_offset']}, "
            info += f"width={sub['width']}, height={sub['height']}, "
            info += f"x_range=({sub['x_offset']}, {sub['x_end'] - 1}), "
            info += f"y_range=({sub['y_offset']}, {sub['y_end'] - 1})\n"
        return info
# # 示例使用
# if __name__ == "__main__":
#     # 定义不同矩阵和划分方式
#     configurations = [
#         {'row': 6144, 'col': 6160, 'subrow': 2, 'subcol': 2},
#         {'row': 9232, 'col': 9216, 'subrow': 2, 'subcol': 8}
#     ]

#     for config in configurations:
#         logger.info(f"Processing matrix {config['row']}x{config['col']} with partition {config['subrow']}x{config['subcol']}")
#         partitioner = MatrixPartitioner(config['row'], config['col'], config['subrow'], config['subcol'])
#         logger.info(partitioner)

#         # 测试一些坐标
#         test_coords = [
#             (0, 0),
#             (3072, 3080),        # 中心区域
#             (6159, 6143),        # 右下角
#             (1152, 4616),        # 边界情况
#             (9215, 9231)         # 右下角最大
#         ]

#         for x, y in test_coords:
#             try:
#                 submatrix = partitioner.find_submatrix(x, y)
#                 range_info = partitioner.find_submatrix_with_range(x, y)
#                 logger.info(f"坐标 (x={x}, y={y}) 位于子矩阵 {submatrix}")
#                 logger.info(f"  详细信息: 子矩阵编号={range_info['submatrix_index']}, "
#                       f"Y范围={range_info['y_range']}, X范围={range_info['x_range']}")
#             except ValueError as e:
#                 logger.info(f"坐标 (x={x}, y={y}) 错误: {e}")
#         logger.info("\n" + "="*80 + "\n")


from astropy.io import fits
from astropy.wcs import WCS
import numpy as np

def ra_dec_to_xy(fits_file, ra, dec):
    """
    将给定的 RA 和 Dec 转换为图像坐标 (x, y)。

    参数:
    fits_file (str): FITS 文件的路径
    ra (float): 右升角（RA），单位为度
    dec (float): 赤纬（Dec），单位为度

    返回:
    tuple: 图像坐标 (x, y)
    """
    # 读取 FITS 文件
    with fits.open(fits_file) as hdul:
        w = WCS(hdul[0].header)  # 获取 WCS 信息
        data = hdul[0].data  # 获取图像数据（可选）

    # 将 RA 和 Dec 转换为像素坐标 (x, y)
    x, y = w.all_world2pix(ra, dec, 0)  # 0 表示使用第一个图层

    return x, y

def ra_dec_to_xy_and_size(fits_file, ra, dec):
    """
    将给定的 RA 和 Dec 转换为图像坐标 (x, y)，并返回 FITS 图像的大小。

    参数:
    fits_file (str): FITS 文件的路径
    ra (str或float): 右升角（RA），可以是 HMS 格式或度（单位为度）
    dec (str或float): 赤纬（Dec），可以是 DMS 格式或度（单位为度）

    返回:
    tuple: 图像坐标 (x, y) 和 FITS 图像的大小 (height, width)
    """
    
    # 将 RA 和 Dec 转换为度（如果是 HMS/DMS 格式）
    if isinstance(ra, str) or isinstance(dec, str):
        coord = SkyCoord(ra=ra, dec=dec, unit=('hourangle', 'deg'), frame='icrs')
        ra = coord.ra.deg
        dec = coord.dec.deg

    # 读取 FITS 文件
    with fits.open(fits_file) as hdul:
        w = WCS(hdul[0].header)  # 获取 WCS 信息
        data = hdul[0].data  # 获取图像数据

        # 获取图像的大小
        height, width = data.shape  # 假设数据是二维的

    # 将 RA 和 Dec 转换为像素坐标 (x, y)
    x, y = w.all_world2pix(ra, dec, 0)  # 0 表示使用第一个图层

    return (x, y), (height, width)